Script for mixed model analysis of Diurnal Saturated Photosynthesis Index DSPI This is a two species analysis and needs to be run in sequence. You can skip the first species but will need to run some of the code located in the Ulva section in order to run linear regressions with growth.

Load the various libraries

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)

Load the dataset

ek_irrad_data <- read.csv("input_data/mean_supersaturation_by_period.csv")

Assign several variables as factors

ek_irrad_data$Run <- as.factor(ek_irrad_data$Run)
ek_irrad_data$Temperature <- as.factor(ek_irrad_data$Temp...C.)
ek_irrad_data$Treatment <- as.factor(as.character(ek_irrad_data$Treatment))
ek_irrad_data$RLC.Order <- as.factor(ek_irrad_data$RLC.Order)

Convert carbon_total from micromol to mol

ek_irrad_data$carbon_mol <- ek_irrad_data$carbon_total / 1e+6

Subset data by species, remove treatment 2.5 for Ulva

ulva <- subset(ek_irrad_data, Species == "ul" & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

ULVA____________________________________________________________

Make a histogram and residual plots of the data for ulva

hist(ulva$carbon_mol, main = paste("Ulva lactuca"), col = "olivedrab3", labels = TRUE)

ulva %>% ggplot(aes(carbon_mol)) +
        geom_histogram(binwidth=0.25, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Run model without interaction between the treatments and temperature - supersat_avg is in minutes Take RLC.Order our of the random effects because causing problems of singularity. R2 is same with or without (+ (1 | RLC.Order))

dspi_model_ulva <- lmer(formula = carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

Plot residuals

hist(resid(dspi_model_ulva))

plot(resid(dspi_model_ulva) ~ fitted(dspi_model_ulva))

qqnorm(resid(dspi_model_ulva))
qqline(resid(dspi_model_ulva))

Check the performance of the model, plot effects, get R2, make table

performance::check_model(dspi_model_ulva)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(dspi_model_ulva)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
##            R2m       R2c
## [1,] 0.3675672 0.8630417
summary(dspi_model_ulva)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: ulva
## 
## REML criterion at convergence: 366.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4455 -0.5733 -0.0052  0.5933  2.8076 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept) 0.075598 0.27495 
##  Run       (Intercept) 0.552863 0.74355 
##  RLC.Order (Intercept) 0.007614 0.08726 
##  Residual              0.175823 0.41931 
## Number of obs: 240, groups:  Plant.ID, 96; Run, 7; RLC.Order, 6
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)    3.394013   0.536212  5.157025   6.330  0.00129 **
## Treatment1     1.374628   0.630693  5.039348   2.180  0.08073 . 
## Treatment2     1.482504   0.630693  5.039348   2.351  0.06511 . 
## Treatment3     1.831997   0.630693  5.039348   2.905  0.03329 * 
## Treatment4     1.881808   0.630693  5.039348   2.984  0.03037 * 
## Temperature27 -0.007210   0.126350  4.925564  -0.057  0.95674   
## Temperature30  0.005544   0.111692 29.335454   0.050  0.96075   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.833                                   
## Treatment2  -0.833  0.991                            
## Treatment3  -0.833  0.991  0.991                     
## Treatment4  -0.833  0.991  0.991  0.991              
## Temperatr27 -0.110 -0.001 -0.001 -0.001 -0.001       
## Temperatr30 -0.105  0.000  0.000  0.000  0.000  0.454
plot(allEffects(dspi_model_ulva))

tab_model(dspi_model_ulva, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  carbon_mol
Predictors Estimates std. Error CI Statistic p df
(Intercept) 3.39 0.54 2.34 – 4.45 6.33 <0.001 229.00
Treatment [1] 1.37 0.63 0.13 – 2.62 2.18 0.030 229.00
Treatment [2] 1.48 0.63 0.24 – 2.73 2.35 0.020 229.00
Treatment [3] 1.83 0.63 0.59 – 3.07 2.90 0.004 229.00
Treatment [4] 1.88 0.63 0.64 – 3.12 2.98 0.003 229.00
Temperature [27] -0.01 0.13 -0.26 – 0.24 -0.06 0.955 229.00
Temperature [30] 0.01 0.11 -0.21 – 0.23 0.05 0.960 229.00
Random Effects
σ2 0.18
Ï„00 Plant.ID 0.08
Ï„00 Run 0.55
Ï„00 RLC.Order 0.01
ICC 0.78
N Run 7
N Plant.ID 96
N RLC.Order 6
Observations 240
Marginal R2 / Conditional R2 0.368 / 0.863

Construct null model to perform likelihood ratio test REML must be FALSE

ulva_dspi_treatment_null <- lmer(formula = carbon_mol ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_dspi_model2 <- lmer(formula = carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_dspi_treatment_null, ulva_dspi_model2)
## Data: ulva
## Models:
## ulva_dspi_treatment_null: carbon_mol ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_dspi_model2: carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                          npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_dspi_treatment_null    7 417.18 441.54 -201.59   403.18          
## ulva_dspi_model2           11 373.07 411.35 -175.53   351.07 52.108  4
##                          Pr(>Chisq)    
## ulva_dspi_treatment_null               
## ulva_dspi_model2          1.309e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_dspi_temperature_null <- lmer(formula = carbon_mol ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_dspi_model3 <- lmer(formula = carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_dspi_temperature_null, ulva_dspi_model3)
## Data: ulva
## Models:
## ulva_dspi_temperature_null: carbon_mol ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_dspi_model3: carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## ulva_dspi_temperature_null    9 369.13 400.45 -175.56   351.13          
## ulva_dspi_model3             11 373.07 411.35 -175.53   351.07 0.0586  2
##                            Pr(>Chisq)
## ulva_dspi_temperature_null           
## ulva_dspi_model3               0.9711

Plots

ulva %>% ggplot(aes(treatment_graph, carbon_mol)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment", y= "DSPI (mol m-2)", title= "E", subtitle = "Ulva lactuca") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", 
                                    "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(0, 10) + stat_mean() + 
        geom_hline(yintercept=5.0, color = "red", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means for DSPI

ulva %>% group_by(Treatment) %>% summarise_at(vars(carbon_mol), list(mean = mean))
## # A tibble: 5 Ă— 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          3.39
## 2 1          4.85
## 3 2          4.96
## 4 3          5.30
## 5 4          5.35
ulva %>% group_by(Run) %>% summarise_at(vars(carbon_mol), list(mean = mean))
## # A tibble: 7 Ă— 2
##   Run    mean
##   <fct> <dbl>
## 1 1      5.91
## 2 2      5.67
## 3 3      5.30
## 4 3.5    4.11
## 5 4      4.18
## 6 6      3.37
## 7 6.5    3.42
ulva %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 7 Ă— 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      643.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 6      626.
## 7 6.5    633.

Add growth rate from other dataset to this one and subset by species

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)

Make a new column for weight change (difference final from initial)

growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100

Subset growth data

gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)

Plot growth vs dspi

ulva_growth_dspi_graph <- ggplot(ulva, aes(x=carbon_mol, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        labs(title = "A", subtitle = "Ulva lactuca", x = "DSPI (mol m-2)", 
             y = "growth rate (%)") + stat_regline_equation(label.x = 5, label.y = 225) + 
        stat_cor(label.x = 5, label.y = 215) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05))
ulva_growth_dspi_graph
## `geom_smooth()` using formula = 'y ~ x'

Plot individual regressions for each treatment vs. growth

ulva_growth_carbon_indiv <- ggplot(ulva, aes(x=carbon_mol, y=growth_rate)) +
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
        geom_smooth(method = "lm", col = "black") + theme_bw() +
        labs(title = "A", subtitle = "Ulva lactuca", x = "DSPI (mol m-2)", 
             y = "growth rate (%)") + 
        stat_cor() +
        #stat_regline_equation() + 
        facet_wrap(vars(Treatment), scales = "free")
ulva_growth_carbon_indiv
## `geom_smooth()` using formula = 'y ~ x'

HYPNEA____________________________________________________________

Subset Hypnea data from larger dataset

hypnea <- subset(ek_irrad_data, Species == "hm" & day1_rlc_time != "11:34:05")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

Make a histogram and residual plots of the data for ulva

hist(hypnea$carbon_mol, main = paste("Hypnea musciformis"), col = "#993333", labels = TRUE)

hypnea %>% ggplot(aes(carbon_mol)) +
        geom_histogram(binwidth=0.25, fill = "#993333", color = "black", size = 0.25, alpha = 0.85) +
        theme_bw()

Run model without interaction between the treatments and temperature - supersat_avg is in minutes Take RLC.Order our of the random effects because causing problems of singularity. R2 is same with or without (+ (1 | RLC.Order))

dspi_model_hypnea <- lmer(formula = carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea)

Plot residuals

hist(resid(dspi_model_hypnea))

plot(resid(dspi_model_hypnea) ~ fitted(dspi_model_hypnea))

qqnorm(resid(dspi_model_hypnea))
qqline(resid(dspi_model_hypnea))

Check the performance of the model, get R2, make table, plot effects

performance::check_model(dspi_model_hypnea)
## Variable `Component` is not in your data frame :/

r.squaredGLMM(dspi_model_hypnea)
##            R2m       R2c
## [1,] 0.1160618 0.8698144
summary(dspi_model_hypnea)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##    Data: hypnea
## 
## REML criterion at convergence: 465.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8895 -0.5737  0.0252  0.4956  4.2986 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Plant.ID (Intercept) 0.09999  0.3162  
##  Run      (Intercept) 0.97668  0.9883  
##  Residual             0.18596  0.4312  
## Number of obs: 286, groups:  Plant.ID, 144; Run, 8
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)     4.16075    0.70521   5.02160   5.900  0.00196 **
## Treatment1      0.29012    0.83424   5.01728   0.348  0.74213   
## Treatment2      0.28747    0.83424   5.01728   0.345  0.74438   
## Treatment2.5    1.22809    1.21529   4.92101   1.011  0.35932   
## Treatment3      0.31804    0.83424   5.01728   0.381  0.71864   
## Treatment4      0.23829    0.83435   5.01995   0.286  0.78660   
## Temperature27  -0.24560    0.09510 115.65047  -2.582  0.01106 * 
## Temperature30  -0.25532    0.09525 116.77735  -2.681  0.00841 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.840                                          
## Treatment2  -0.840  0.994                                   
## Treatmnt2.5 -0.577  0.488  0.488                            
## Treatment3  -0.840  0.994  0.994  0.488                     
## Treatment4  -0.840  0.994  0.994  0.487  0.994              
## Temperatr27 -0.067  0.000  0.000  0.000  0.000  0.000       
## Temperatr30 -0.067  0.000  0.000  0.000  0.000  0.001  0.499
plot(allEffects(dspi_model_hypnea))

tab_model(dspi_model_hypnea, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
  carbon_mol
Predictors Estimates std. Error CI Statistic p df
(Intercept) 4.16 0.71 2.77 – 5.55 5.90 <0.001 275.00
Treatment [1] 0.29 0.83 -1.35 – 1.93 0.35 0.728 275.00
Treatment [2] 0.29 0.83 -1.35 – 1.93 0.34 0.731 275.00
Treatment [2.5] 1.23 1.22 -1.16 – 3.62 1.01 0.313 275.00
Treatment [3] 0.32 0.83 -1.32 – 1.96 0.38 0.703 275.00
Treatment [4] 0.24 0.83 -1.40 – 1.88 0.29 0.775 275.00
Temperature [27] -0.25 0.10 -0.43 – -0.06 -2.58 0.010 275.00
Temperature [30] -0.26 0.10 -0.44 – -0.07 -2.68 0.008 275.00
Random Effects
σ2 0.19
Ï„00 Plant.ID 0.10
Ï„00 Run 0.98
ICC 0.85
N Run 8
N Plant.ID 144
Observations 286
Marginal R2 / Conditional R2 0.116 / 0.870

Construct null model to perform likelihood ratio test REML must be FALSE

hypnea_dspi_treatment_null <- lmer(formula = carbon_mol ~ Temperature + (1 | Run) + (1 | Plant.ID) , data = hypnea, REML = FALSE)
hypnea_dspi_model2 <- lmer(formula = carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea, REML = FALSE)
anova(hypnea_dspi_treatment_null, hypnea_dspi_model2)
## Data: hypnea
## Models:
## hypnea_dspi_treatment_null: carbon_mol ~ Temperature + (1 | Run) + (1 | Plant.ID)
## hypnea_dspi_model2: carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_dspi_treatment_null    6 466.64 488.57 -227.32   454.64          
## hypnea_dspi_model2           11 474.26 514.48 -226.13   452.26 2.3719  5
##                            Pr(>Chisq)
## hypnea_dspi_treatment_null           
## hypnea_dspi_model2             0.7956
hypnea_dspi_temperature_null <- lmer(formula = carbon_mol ~ Treatment + (1 | Run) + (1 | Plant.ID), data = hypnea, REML = FALSE)
hypnea_dspi_model3 <- lmer(formula = carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID), data = hypnea, REML = FALSE)
anova(hypnea_dspi_temperature_null, hypnea_dspi_model3)
## Data: hypnea
## Models:
## hypnea_dspi_temperature_null: carbon_mol ~ Treatment + (1 | Run) + (1 | Plant.ID)
## hypnea_dspi_model3: carbon_mol ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID)
##                              npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_dspi_temperature_null    9 479.49 512.40 -230.75   461.49          
## hypnea_dspi_model3             11 474.26 514.48 -226.13   452.26 9.2292  2
##                              Pr(>Chisq)   
## hypnea_dspi_temperature_null              
## hypnea_dspi_model3             0.009906 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plots

hypnea %>% ggplot(aes(treatment_graph, carbon_mol)) + 
        geom_boxplot(size=0.5) + 
        geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
        labs(x="Treatment", y= "DSPI (mol m-2)", title= "F", subtitle = "Hypnea musciformis") + 
        scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", 
                                    "18ppt/53umolN", "11ppt/80umolN")) + 
        ylim(0, 10) + stat_mean() + 
        geom_hline(yintercept=5.0, color = "red", size = 0.5, alpha = 0.5) +
        scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
        theme_bw() +
        theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
              plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

Summarize the means for DSPI

hypnea %>% group_by(Treatment) %>% summarise_at(vars(carbon_mol), list(mean = mean))
## # A tibble: 6 Ă— 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          3.99
## 2 1          4.45
## 3 2          4.45
## 4 2.5        5.22
## 5 3          4.48
## 6 4          4.41
hypnea %>% group_by(Run) %>% summarise_at(vars(carbon_mol), list(mean = mean))
## # A tibble: 8 Ă— 2
##   Run    mean
##   <fct> <dbl>
## 1 1      5.55
## 2 2      5.29
## 3 3      4.10
## 4 3.5    3.07
## 5 4      3.37
## 6 7      5.22
## 7 8      4.12
## 8 9      3.87
hypnea %>% group_by(Run) %>% summarise_at(vars(day_length_avg), list(mean = mean))
## # A tibble: 8 Ă— 2
##   Run    mean
##   <fct> <dbl>
## 1 1      654.
## 2 2      645.
## 3 3      626.
## 4 3.5    621.
## 5 4      611.
## 6 7      688.
## 7 8      635.
## 8 9      636.

For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/29/21 because it was white and also looked dead

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

Plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_carbon_graph <- ggplot(hypnea, aes(x=carbon_mol, y=growth_rate)) + 
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
        geom_smooth(method = "lm", col = "black") + theme_bw() + 
        xlim(2, 8) + ylim(-100, 190) + 
        labs(title = "B", subtitle = "Hypnea musciformis", x = "DSPI (mol m-2)", 
             y = "growth rate (%)") + 
        stat_regline_equation(aes(size = 14), label.x = 6, label.y = -75) + 
        scale_color_manual(values = c("#57330F", "#22773E", "#9970C2", "#FF3333", "#FF9933", "#99FF33")) +
        stat_cor(aes(size = 14), label.x = 6, label.y = -85) +
        theme(plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05, size = 14), 
              plot.subtitle = element_text(face = "italic", vjust = -20, hjust = 0.05, size = 14))
hypnea_growth_carbon_graph
## `geom_smooth()` using formula = 'y ~ x'

Plot individual regressions for each treatment vs. growth

hypnea_growth_carbon_indiv <- ggplot(hypnea, aes(x=carbon_mol, y=growth_rate)) +
        geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
        geom_smooth(method = "lm", col = "black") + theme_bw() +
        labs(title = "B", subtitle = "Hypnea musciformis", x = "DSPI (mol m-2)", 
             y = "growth rate (%)") + 
        stat_cor() +
        #stat_regline_equation() + 
        facet_wrap(vars(Treatment), scales = "free")
hypnea_growth_carbon_indiv
## `geom_smooth()` using formula = 'y ~ x'